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Characterizing structural inhomogeneity is an essential step in understanding the mechanical 
response of amorphous materials. We introduce a threshold-free measure based on the field of 
vectors pointing from the center of each particle to the centroid of the Voronoi cell in which the 
particle resides. These vectors tend to point in toward regions of high free volume and away from 
regions of low free volume, reminiscent of sinks and sources in a vector held. We compute the local 
divergence of these vectors, where positive values correspond to overpaeked regions and negative 
values identify underpaeked regions within the material. Distributions of this divergence are nearly 
Gaussian with zero mean, allowing for structural characterization using only the moments of the 
distribution. We explore how the standard deviation and skewness vary with packing fraction for 
simulations of bidisperse systems and hnd a kink in these moments that coincides with the jamming 
transition. 
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In glassy liquids and disordered solids, heterogeneities 
in local structure correlate with heterogeneous particle 
rearrangement dynamics arising from thermal fluctua¬ 
tions or applied mechanical load HHS]. Gharacterizing 
local structural heterogeneity is therefore important in 
experiment, for example via the network of contacts and 
force chains IMO], and as a step in understanding ther¬ 
mal and mechanical response. A simple and physically- 
appealing measure of local structure that forms the basis 
of historically-important theories of glassy dynamics and 
plasticity is free volume Regions that are un¬ 

derpacked have a larger local free volume, and therefore 
ought to rearrange or yield more easily. Though intu¬ 
itive, free volume is inherently a concept based on hard 
spheres and only applies at densities below jamming. 

Here we introduce a generalization of the concept of 
free volume that derives from the radical Voronoi net¬ 
work and hence applies in a consistent way to particles 
interacting via any inter-particle potential at any density. 
Our measure, Q/e, is inspired by the observation that the 
center of a particle center deviates from the centroid of 
the corresponding Voronoi cell when the configuration is 
disordered. In two dimensions, it is defined as 

Qk = {V ■c){Ak/A), (1) 

where c is the interpolated field of vectors that point 
from particle centers to the corresponding Voronoi cell 
centroids, the divergence is taken over a Delaunay trian¬ 
gle k with area and A is the average of all Aj^ within 
the a packing. By construction, Qk is dimensionless and 
has zero mean. It is sensitive to local structural hetero¬ 
geneity and - though purely geometrical - has a clear 
physical interpretation: positive/negative values respec¬ 
tively correspond to overpacked/underpacked regions. In 
addition to establishing a statistical correlation between 
Qk and local relative free volume, we find that the dis¬ 
tribution of Qk values over a packing is nearly Gaussian, 


with mode and median nearly equal to the mean (zero); 
hence it may be well-described by just the standard de¬ 
viation and the skewness. 

As an illustration, we calculate Qk for a system of soft 
disks at a series of packing fractions that are widely var¬ 
ied, from the dilute limit to well above the jamming tran¬ 
sition. Morse and Gorwin m have recently identified 
geometrical features similarly based on Voronoi tessella¬ 
tions that exhibit singularities at the jamming transition. 
Here we find the standard deviation and skewness of the 
Q/c-distribution also exhibit kinks at the transition. Thus 
there is a signature of the jamming transition in Q/., a 
geometrical quantity with clear physical relevance. 



FIG. 1: (color online) Particle packing from (a) simula¬ 
tion and (b) experiment [161 (13? with superimposed rad¬ 
ical Voronoi tesselation (blue) and Delaunay triangulation 
(green). Also shown are vectors Cp (magenta) that point from 
each particle center (red dot) to the centroid of its Voronoi 
cell. In (b), these are elongated by a factor of 8 for each of 
visualization. 


We begin by using the software package voro++ [18] to 
determine the radical Voronoi tessellation, a space-filling 
generalization of the standard Voronoi construction to 
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polydisperse systems. In this framework, each cell edge 
is determined from two adjacent circles, and is given by 
the locus of points from which tangent lines drawn to 
both circles have the same length. If all particle radii 
are equal, the standard Voronoi tessellation is recovered. 
Fig.[T] shows an example, where the Voronoi tessellation 
is overlaid on an image of particles. Any two particles 
with a shared Voronoi cell edge are defined as neighbors, 
and from this, we generate a generalized Delaunay trian¬ 
gulation by connecting groups of three mutual neighbors 
into triangles, as shown in green in Fig. 

The position of a particle within its Voronoi cell is an 
indicator of local variation in the packing. This moti¬ 
vates consideration of a local anisotropy vector, C^, for 
each particle p, that points from the center of the par¬ 
ticle to the centroid of its Voronoi cell, as shown by the 
magenta arrows in Fig. In monodisperse crystalline 
packings with a single particle per unit cell, each par¬ 
ticle center and the corresponding Voronoi cell centroid 
coincide; therefore = 0 for all p, consistent with the 
idea that is a measure of the structural anisotropy. 
This vector is one of several Minkowski functionals m 
associated with a Voronoi cell, many of which have been 
used to describe packing heterogeneity [201 EH. Addi¬ 
tionally, it has been discussed in the context of structure 
in liquids [22l [23], and has been found to be correlated 
with particle motion [221121]. As might be expected, 
points in the direction of excess free volume, indicating 
the direction in which the particle has the most space 
to move. However, local spatial variations of this vector 
have not been previously explored. 

Fig. Et shows a typical example of vectors for sev¬ 
eral particles in a bidisperse packing with particle size 
ratio of 3:4 and hard-sphere interactions. Vectors tend 
to point in toward locally less well-packed and away from 
locally more well-packed regions of the packing, reminis¬ 
cent of sinks and sources in a vector field. Therefore it is 
natural to consider the divergence of a field defined by in¬ 
terpolating the Cp vectors over a local region. We choose 
Delaunay triangles as the local regions over which to per¬ 
form interpolation and differentiation of the vectors. 
This allows us to use the framework of the constant strain 
triangle of finite element analysis [25] to find local spatial 
variations of the Cp vectors. In particular, each triangle 
is treated independently and we assume that the asso¬ 
ciated Cp vectors define a vector field c = (cx^Cy) that 
varies linearly over each triangle: 

Cx{x, y) =dx^ dxxX + dxyV, 

^y{Xi - dy “h dyxX “h dyyy. (2) 


vector (dx^dy) is equivalent to (cx^Cy)^ where the av¬ 
erages are taken over the triangle vertices. The tensor 
dij = dcijdxj is independent of the origin location and 
describes the spatial variation of the c field over the tri¬ 
angle, with the divergence given by Ti{dij). 



Q, = (V.c)(A,/A) 
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FIG. 2: (color online) (a) Probability density of normalized 
divergence of center-to-centroid vectors for the experimen¬ 
tal packing of bidisperse hard spheres shown in Fig. (cir¬ 
cles) and best-fit Gaussian (solid curve). Qk > 0 regions are 
more tightly packed than their surroundings, hence we call 
these regions overpacked. Qk < 0 regions are more loosely 
packed than their surroundings, and are therefore labeled un¬ 
derpacked. The inset shows the probability density and Gaus¬ 
sian fit with a logarithmic y-axis, highlighting the deviation of 
the data from Gaussian, particularly in the negative tail, b) 
PDFs of triangle-based free area fraction, 1 —and Voronoi- 
based free area fraction, 1 — 0p, for the same packing. These 
distributions are more complicated and less intuitive than for 
Qfc. Inset: Qk correlates well with a fractional deviation from 
the global packing fraction, and therefore is similar to a rel¬ 
ative free area. All dashed lines indicate the median of the 
data set of the same color. 


For each triangle the six constants dx^dy^dxx^dxy^dyx^ 
and dyy can be determined by evaluating Eq. Q at 
the triangle vertices and inverting the resulting matrix 
equation. If the triangle coordinates are shifted so that 
the centroid of the triangle is located at the origin, the 


From the divergence theorem, all contributions from 
interior particles cancel upon performing the sum Qk 
over all triangles in a packing, leaving only contributions 
from the boundary particles. This results in {Qk) = 0 for 
infinite systems and for systems with periodic boundary 
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conditions. The interpretation of Qk is then physically 
intuitive: triangles with Qj^ < 0 are less well packed than 
their surroundings, so we label these regions as under- 
packed^ while triangles with > 0 are locally more 
tightly packed than their surroundings, so we refer to 
these regions as overpacked. Thus Qk is a measure of 
relative free volume, and, as shown below, the statistical 
correlation is quantitative. 

Fig. shows the probability density of Qk for the ex¬ 
perimental bidisperse system shown in Fig.[^. The prob¬ 
ability density is nearly Gaussian with a small, positive 
mean resulting entirely from the finite boundaries of the 
packing. Thus we may characterize packings to a great 
degree just from the standard deviation and skewness of 
the Qk distribution. To the extent that the probability 
density is truly Gaussian and the Qk values are spatially 
uncorrelated [26] , the packing is random in a very simple 
sense. But since adjacent triangles share two vectors, 
there must be at least short-range correlations. Never¬ 
theless, Qk is closer to a Gaussian random variable than 
any other structural quantity previously used to charac¬ 
terize random packings. Furthermore, deviations from 
Gaussianity (e.g. underpacked particles in the tail of the 
distribution) are likely to have important physical conse¬ 
quences [27] . 

To build intuition, we now compare Qk with the lo¬ 
cal area fraction. Fig. shows the probability density 
for two standard measures, based on particle area per 
Voronoi cell and per Delaunay triangle. The probability 
densities for these quantities have irregular complicated 
shapes, where the median differs significantly from the 
mode for each distribution. There is no clear feature 
demarcating under- versus over-packed regions. Never¬ 
theless, the magnitude of Qk is on the order of - and 
statistically correlated with - the relative free area de¬ 
fined as {(j)k — where <pk is the triangle-based area 
fraction and 0 is the global area fraction of the sample. 
This is demonstrated by the contour plot of relative free 
area versus Qk in the inset of Fig.[^. Similar plots in [26] 
show that good correlation holds at all packing fractions, 
but less so for dilute systems due to the development of 
long tails away from the heart of the distributions. We 
may thus consider the size of Q/c as a semi-quantitative 
indication of local free volume relative to the average 
packing. 

For the remainder of the paper, we use the Qk distri¬ 
bution as a tool to characterize structure versus packing 
fraction for simulated systems in two-dimensions. Static 
packings are created using four different protocols. For 
the first, a large number N = 5000 to = 80,000 of 
points are placed at random in a box. For each we 
find that the average moments of Qk over 200 configu¬ 
rations are independent of N. For the second protocol, 
we generate several packings of non-overlapping monodis- 
perse disks. Here, each disk has a radius equal to 1, and a 
proposed new disk is only accepted and placed in the box 


if it does not overlap with any existing disks. These pack¬ 
ings contain anywhere from N = 1000 to N = 160,000 
disks, spanning area fractions from 0 0.003 to 0 0.5. 

At least 40 packings are generated for each value of (j). 

For the final two protocols, we numerically generate 
systems composed of = 2048 soft repulsive disks with 
mass M. Two particles i and j interact with the pairwise 
potential 
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where rjj is the distance between the particle centers, Rj 
and Rj are the particle radii, Q{x) is the Heaviside step 
function, and e sets the energy scale. To prevent crystal¬ 
lization, we use a 50 — 50 mixture of particles with a size 
ratio of 1.4. The disks are initially placed at random in 
a periodic simulation box with a packing fraction 0, and 
are then allowed to move according to two different pro¬ 
cedures. The units of length, mass, and energy are 2i?avg5 
M, and e respectively, where R^vg = 
average particle radius. 

For the third protocol, thermalized configurations are 
generated at a very low temperature using molecular dy¬ 
namics simulations at constant NVT, performed using 
LAMMPS [28|. Beginning at a temperature of Tgtart = 
0.05, we slowly cool the system to T = 10“^ over 5 x 10^ 
time steps. The system then remains at T = 10“^ for 
an additional 10^ time steps. The fourth and final proto¬ 
col corresponds to an infinitely fast quench from infinite 
to zero temperature. Beginning with the initial random 
configuration, we minimize the total energy to a local 
minimum using the FIRE algorithm [29]. Each protocol 
was repeated 500 times at each packing fraction. 

Eor all final configurations, the Qk distributions and 
moments are computed. The low-0 behavior of the stan¬ 
dard deviation is emphasized in the main plot of Eig.[^. 
The results for protocols 2 — 4 appear to converge nicely 
to the value 0.5746 ± 0.0002 found from the random 
point patterns of protocol 1. Eig. shows a similar 
convergence of the skewness consistent with the value 
—2.086 ± 0.0002 obtained for point patterns. This vali¬ 
dates the protocol methods, and serves to establish the 
low-0 “ideal gas” limiting behavior of the Qk distribu¬ 
tions. 

As 0 increases away from zero, the moments of the Qk 
distribution change in a protocol-dependent fashion. As 
seen, the thermalized configurations are closer than the 
rapidly-quenched configurations to the randomly-placed 
non-overlapping sphere configurations. The Qk distribu¬ 
tions all become narrower with increasing 0, as shown 
by the standard deviation plot. And as judged from the 
skewness plot, the distributions generally become more 
Gaussian - as mentioned earlier. However the quenched 
configurations show an initial increase in the skewness 
magnitude before decreasing towards zero. 
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FIG. 3: (color online) Standard deviation (a) and skewness 
(b) of the Qk distributions versus global packing fraction 0. 
Each data point represents the average value of the speci- 
hed moment over several conhgurations; the associated un¬ 
certainty is smaller than the symbol. The main plots em¬ 
phasize the low 0 behavior, which approaches the same con¬ 
stant value for all preparation protocols shown: inhnitely-fast 
quench (x), thermal (o), randomly placed non-overlapping 
monodisperse circles (+), and random point patterns with 
0 = 0 (dashed line). Insets show behavior near 0c- Rescaled 
distributions of 0c are plotted as solid curves, with color indi¬ 
cating preparation protocol from which they were determined. 
The moments of Qk have a kink that coincides with the peak 
in the respective 0c distribution. 


The insets in Fig. show zoom-ins of the moments 
near the critical packing fraction 0c, where the sys¬ 
tems become jammed. For our simulations, there is 
a distribution of 0c values due to the finite system 
size [30l [3T] . Each protocol produces a different distribu¬ 
tion [32l[33]. Here, the average and standard deviation 
are 0c = 0.8409 ±0.0012 for the quenched protocol and 
0c = 0.8465 ± 0.0005 for the thermalized protocol. As 
expected, the 0c values are smaller and more widely dis¬ 
tributed for quenched configurations, since thermalized 
configurations have more opportunity to relax [32]. The 
unnormalized 0c distributions are individually rescaled 
to reach the respective data curve in the insets of Fig. [^ 
in order to mark the jamming transitions. 


The key striking result evident in the Fig. insets is 
that a signature of the jamming transition exists in the 
Qk distributions. Namely, the standard deviation and 
the skewness both show a kink where the (pc distribu¬ 
tions are peaked. For the quenched protocol, the skew¬ 
ness kink is smallest and may deviate from slightly 0c. 
For the thermalized protocol, the behavior near 0c is con¬ 
siderably more dramatic. In particular, the kinks are ex¬ 
tremely pronounced in that the derivatives of standard 
deviation and of skewness versus 0 actually change sign 
on opposite sides of the transition. Furthermore, there 
is non-monotonic behavior below the transition, with the 
standard deviation and the skewness exhibiting a mini¬ 
mum and maximum, respectively, at a packing fraction 
a few percent below 0c. As a measure of static struc¬ 
ture versus 0, the Qk distribution is thus even capable 
of signaling a precursor that the onset of jamming is im¬ 
minent. Since this happens for the thermalized but not 
the quenched configurations, it could be related to the 
dynamical hard-sphere glass transition; however, the ex¬ 
trema are at different 0 values. Note that the differences 
in the trends for the two protocols implies that details of 
the local structure are sensitive to protocol even though 
other quantities that are singular at the jamming transi¬ 
tion, such as the average contact number, scale the same 
way with increasing pressure for packings prepared using 
different protocols [32]. 


In conclusion, we have shown that Qk is a geometrical 
measure of structure, with the physical meaning of a rela¬ 
tive free volume, which displays a strong signature of the 
jamming transition. It is particularly easy to interpret 
Qk in terms of overpacked and underpacked regions since 
the Q/c-distribution has zero mean, by construction, and 
is nearly Gaussian for non-dilute systems. Though all our 
examples are two-dimensional, the concept of Qk may 
be extended to any dimension by appropriate Voronoi 
construction. For thermal and sheared systems, there 
is longstanding interest in identifying structural features 
that lead to dynamical activity such as heterogeneous 
particle rearrangements and shear bands. The correla¬ 
tion of Qk with dynamics, as well as with structural pre¬ 
dictors of rearrangements found by machine learning [6] , 
may now be studied. This could give new meaning to 
the concept of “free volume,” which has been assumed 
in many theories to affect dynamics and thereby control 
the glass transition and glassy rheology [T3] . 
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Supplemental figures for “Divergence of Voronoi cell anisotropy vector: A 
threshold-free characterization of local structure in amorphous materials” 
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FIG. Al: (a) Center-to-centroid vectors C for an experimental packing, like the one shown in Fig. 1 of the main text, elongated 
by a factor of 8 for ease of visualization, (b) Resulting Qk field for the same region of the same packing. More loosely packed 
regions appear green, while locally more well packed regions are blue. 
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FIG. A2: Correlation between the local relative free area and the area-weighted divergence of the center-to-centroid vector field, 
for numerous global area fractions. Here, the packings are numerically generated by the thermalization protocol at T = 10~^. 
The blue line is the same in all plots, y — 2x. The right column shows a range of packing fractions spanning the jamming 
transition at (j)c = 0.8465 ± 0.0005. Throughout this range the relative free area is highly correlated with Q^, and is twice as 
large on average. The left column shows a broader range of cj). Far below (/)c, the relative free area is still roughly 2Qk on 
average, but the correlation in the tails becomes progressively weaker. 
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FIG. A3: Correlation between the local relative free area and the area-weighted divergence of the center-to-centroid vector 
field, for numerous global area fractions. Here, the packings are simulated by rapid quench from T = oo to T = 0, and the 
jamming transition is at (j^c — 0.8409 zb 0.0012. As seen in the previous plot, the relative free area is about 2Qk on average, at 
all 0, but with tails that weaken the correlation upon dilution. 

























9 


Qfe = (V • c){Ak/A) 


-1-1-1-1-r 


—1-1-1-1— 

Quenched 

XXxXxx XX w 

Thermalized 



^^^^ 

Ai 




Dk = \/ -c 





(/) 

C/) 

0 

C 


0 


w 



0.80 0.82 0.84 0.86 0.88 

0 



0.80 0.82 0.84 0.86 0.88 


0 


FIG. A4: Moments for both area-weighted (Qfc, left column) and non-area-wieghted {Dk, right column) divergence of the 
interpolated center-to-centroid vectors. The modes for both Qk and Dk are comparable, and there is no noticeable feature in 
either curve that coincides with the jamming transition. The mean of Qk is zero by construction, while the mean of Dk is 
non-zero. There is a small feature in the mean of Dk that coincides with each transition. The standard deviations are very 
similar, and both have a kink at the transition. Both skewnesses have kinks at the transition, though the skewnesses are more 
negative and the kinks are more pronounced in Qk. 
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Qfc = (V • c){Ak/A) 


FIG. A5: Probability density functions for experimental (black circles), quenched (green), and thermalized (blue) data sets, 
all of which have packing fraction, 0 0.81. The linear y-axis (a) highlights differences in the peaks of the distributions, 

and the logarithmic y-axis (b) emphasizes the differences in the tails of the distributions. The experimental data more closely 
resembles the quenched data, though both numerically-generated packings have sharper peaks and narrower tails than the 
experiment. Some of these differences are due to preparation protocol, and some may be due to differing particle size ratios 
(3:4 in experiment, 1:1.4 in simulation). 
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FIG. A6: Normalized spatial autocorrelation of four different measures of local free area, at equal times, versus radial distance 
divided by average particle radius. These are determined by linearly interpolating each structural measure at every point in 
time, computing the equal-time two-dimensional autocorrelation of the mean-subtracted interpolations, and normalizing by the 
value for zero shift. The result is then averaged over angular shifts at each time, and finally averaged over all times. In all 
cases, correlations are fairly short-ranged, vanishing by Ar ~ 5 particle radii. The negative correlation in both Dk and Qk is 
expected as both are measures of free area relative to the immediate surroundings. 
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FIG. A7: Some methods of structural characterization start with a contact or force network. However, this information is not 
always precisely known, as is the case in for the experimental system shown here. Each line is drawn from one particle center to 
an adjacent particle, and the grayscale shading of these links indicates the distance between particle surfaces. In experimental 
systems such as this one, there is uncertainty in precise particle locations and size, requiring the choice of a separation threshold 
to decide which particles are actually in contact. Lighter colors represent larger separations, so these connections are the first to 
be pruned if the separation threshold is decreased. Very small changes in this threshold result in large changes in the resulting 
network, making this approach unreliable if contact information is imperfectly known as in experiment. Voronoi-based measures 
are much more robust to small uncertainties in particle sizes and locations, making them useful in a wider variety of systems. 




























